Simulation of mode II unconstrained fracture path formation coupled with continuum anisotropic damage propagation in shale

The objective of this work is to simulate mode II multi-scale fracture propagation in shale by coupling a continuum anisotropic damage model with a Cohesive Zone Model (CZM). The Continuum Damage Mechanics – based Differential Stress Induced Damage (DSID) model is used to represent micro-scale crack evolution. DSID parameters were calibrated against pre-peak points of stress/strain curves obtained experimentally during triaxial compression tests performed in Bakken shale. A bilinear CZM is employed to represent macroscale fracture propagation. We calculated the effective shear modulus of a continuum that contained a distribution of parallel cracks according to the DSID model (which does not account for crack interactions) and according to Kachanov’s micromechanical model (which accounts for crack interactions). Simulations confirmed that above a crack density or damage of 0.3, crack interactions could not be ignored, and we used that threshold to define the transition between continuum damage propagation and discrete fracture propagation and subsequently, to calibrate the shear cohesive strength of the CZM. The CZM cohesive energy release rate was determined by calibrating a numerical model of triaxial compression test against experimental data obtained on Bakken shale. The cylindrical sample was modeled with a CZM to pre-define an inclined cohesive fracture, and the DSID model was assigned to the surrounding elements. We used our calibrated CZM-DSID model to simulate a biaxial compression test in plane strain. Results clearly show that the proposed modeling strategy not only allows simulating the advancement of macro-fracture tips, but also captures the inception and growth of micro-cracks that form damaged zones, as well as the transition between smeared damage and discrete fracture. 2. ANISOTROPIC DAMAGE MODEL 2.1. Theoretical Framework of the DISD model We study the microscale crack propagation with the DSID model (Xu and Arson, 2015). The free enthalpy is expressed in terms of elastic energy and damage induced additional energy, in which the damage variable is defined as a second order tensor Ω to represent crack evolution inside an REV. Gs(σ ,Ω) = 1 2 σ : S0 :σ + a1TrΩ Trσ ( )2 + a2Tr σ ⋅σ ⋅Ω ( ) +a3TrσTr Ω⋅σ ( ) + a4TrΩTr σ ⋅σ ( ) (1) Where ai are material parameters. S0 is the reference (undamaged) compliance tensor. The thermodynamic conjugated relationships are the following: ε E = ε el + ε ed = ε − ε id = ∂Gs ∂σ = 1+ν0 E0 σ − ν0 E0 Trσ ( )δ + 2a1 TrΩTrσ ( )δ + a2(σ ⋅Ω +Ω⋅σ )+)a3 Tr Ω⋅σ ( )δ + Trσ ( )Ω ⎡⎣ ⎤⎦ + 2a4 TrΩ ( )σ Y = ∂Gs ∂Ω = a1 Trσ ( )2 δ + a2σ ⋅σ +a3 Trσ ( )σ + a4Tr σ ⋅σ ( )δ (2) In which 0 E and 0 ν are Young’s modulus and Poisson ratio of initial undamaged material. As shown in Fig 1, ε el is the purely elastic strain, ε ed is the elastic damageinduced strain that result from the degradation of mechanical stiffness, and id ε is the irreversible strain. The damage criterion is a Drucker Prager yield function expressed in terms of energy release rate instead of stress, which allows predicting the evolution of damage with deviatoric stress: fd = J * −α I * − k J * = 1 2 P1 : Y1 3 I δ ⎛ ⎝⎜ ⎞ ⎠⎟ : P1 : Y1 3 I δ ⎛ ⎝⎜ ⎞ ⎠⎟ I = P1 : Y ( ) :δ P1 σ ( )= H σ p ( )− H −σ p ( ) ⎣ ⎦ p=1 3 ∑ n p ⊗ n p ⊗ n p ⊗ n p k = C0 −C1Tr(Ω) (3) Where C0 is the initial damage threshold, 1 C is an isotropic hardening variable, ( ) H ⋅ is the Heaviside function, and σ p is the p-th principal stress. In order to satisfy Clausius-Duhem inequality, the damage potential is defined as: gd = 1 2 P2 : Y ( ) : P2 : Y ( ) P2 σ ( )= H maxq=1 3 σ p ( )−σ p ⎣ ⎦ p=1 3 ∑ n p ⊗ n p ⊗ n p ⊗ n p (4) The flow rule is associated for the irreversible damage strain, and non-associated for damage, as follows: ! ε id = ! λd ∂ fd ∂σ = ! λd ∂ fd ∂Y : ∂Y ∂σ ! Ω = ! λd ∂gd ∂Y (5) Fig. 1. Energy components in the DSID model Table 1. Calibrated DSID parameters Parameters Units Value